Solid-solid phase transition in hard ellipsoids 
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Abstract 

We present a computer simulation study of the crystalline phases of hard ellipsoids of revolution. 
A previous study [Phys. Rev. E, 75, 020402 (2007)] showed that for aspect ratios a/b > 3 the 
previously suggested stretched-fcc phase [Mol. Phys., 55, 1171 (1985)] is unstable with respect 
to a simple monoclinic phase with two ellipsoids of different orientations per unit cell (SM2). In 
order to study the stability of these crystalline phases at different aspect ratios and as a function 
of density we have calculated their free energies by thermodynamic integration. The integration 
path was sampled by an expanded ensemble method in which the weights were adjusted by the 
Wang-Landau algorithm. We show that for aspect ratios a/b > 2.0 the SM2 structure is more 
stable than the stretched-fcc structure for all densities above solid-nematic coexistence. Between 
a/b = 1.55 and a/b = 2.0 our calculations reveal a solid-solid phase transition. 
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I. INTRODUCTION 



Suspensions of hard particles (i.e. particles that interact via an infinitely strong, repulsive 
excluded-volume interaction potential) have been successfully used as model systems for the 
statistical mechanics of liquids and solids for more than half a century. For this class of 
system phase transitions are entropy rather than enthalpy driven, and the relevant control 
parameters are the particle shape and concentration rather than temperature. Hard ellip- 
soids are a simple model for systems whose macroscopic properties depend on the interplay 
of positional and orientational entropy such as liquid crystals Q, [3, 3, li] and orientational 
glasses [5|, la, [Z| 

In recent years it has been shown by computer simulations and experiments that ran- 



domly packed arrangements o 



close packing of spheres 



rard ellipsoids can reach densities much higher than random 



At certain aspect ratios, random packing of ellipsoids 
can even reach densities almost as high as the closest crystalline packing of spheres js]. 
However, this does not imply that random packing of ellipsoids is as dense as their densest 
known crystalline packing. In 2004, Donev and co-workers introduced a family of crystalline 
packings of ellipsoids, which reach a packing fraction of rj ~ 0.7707 ll|] (as compared to 
7] = ~ 0.7405 for the fee packing of spheres and stacking variants thereof). 

Inspired by this study, we re-examined the phase diagram of hard ellipsoids |12J . We found 
that the stretched fec-phase, which had before been assumed to be the stable crystalline 



phase 



131 ]. was unstable with respect to a different crystalline phase. The more stable 



structure has a simple monoclinic unit cell containing two ellipsoids of unequal orientation 
(SM2)(cf. Fig [I]). The packings constructed by Donev and co-workers 11] are a special case 
of SM2 (the infinite-pressure limit). 

At that time we did, however, not compute free energy differences between SM2 and 
stretched-fec. In the present article we report on Monte Carlo simulations in which SM2 and 
stretched-fee are connected to their respective harmonic crystals ("Einstein crystals") via 
thermodynamic integration, and hence their free energies are determined. In order to sample 
the thermodynamic integration pathway, we adapted the Wang-Landau algorithm In 
the original Wang-Landau scheme a flat histogram of the internal energy is constructed. Here 
we constructed a flat histogram of the coupling parameter that couples the hard ellipsoid 
model to the Einstein crystal, instead. 
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Figure 1: Unit cell of SM2 [12( with a/b = 3 (color online). The open circles indicate the centers 
of the two ellipsoids which form the basis. The yellow (light gray) ellipsoid is at the origin, the 
green (dark gray) one is at |(a + b). The cell is monoclinic. f3 is the soft degree of freedom. 
Part c) shows the cell at close packing (the infinite-pressure limit), where it is an instance of the 



family of packings introduced by Donev et al. [ll|]. Note the indicated right angle and the resulting 
symmetry about the 6c-plane in this case. 

II. METHOD OF COMPUTATION 

In order to determine which of two phases is thermodynamically more stable, one com- 
pares their relevant thermodynamic potentials, e.g. in the case of constant particle number 
N, volume V and temperature T their free energies F. Within a MC simulation, however, 
for most models it is impossible to compute F because of its direct connection to the acces- 
sible phase space volume (q N ,p N )- To solve this problem the method of Thermodynamic 
Integration (TI) flol . fisj ] is commonly used, in which the free energy difference between 
the system of interest and a reference system can be calculated by introducing an artificial 
external potential U, such that 

Here, the parameter ( e [0, 1] links the interaction potential of the system of interest U Bys = 
U(( = 0) to the potential of the reference system C/ ref = U(( = 1) by 

U(q N ; C) = (1 - C)U sys (q N ) + (U ie{ (q N ) . (2) 

During a typical integration, U SJS is gradually switched on and at the same time U re { is 
gradually switched off. However, the hard-core interaction of the ellipsoids does not allow 
for a gradual change. Therefore, U sys is imposed in a first step, and then LV re f is gradually 



switched off in a second. With this procedure, the free energy of the system can be calculated 
as F sys = F re f + AFi + AF 2 , where the subscripts refer to the two steps just described. 

AF 1 = -ln(exp[-/3[/ sys ]) c=1 , (3) 

where U sys (q N ) is here the overlap potential of the ellipsoids, and the configuration q N 
consists of positions and orientations q N = {r N , 8 N }. (. . .)^ refers to the ensemble average 
where the potential is parameterized by Eq. [21 AF 2 will be discussed together with the free 
energy of the reference state in the following paragraph. 

From Eq. [I] it is obvious that F re f needs to be known from other sources, e.g. by analytical 
computation and that no phase transition may occur during the integration process. In order 
to construct such a reference system, we consider a system of hard ellipsoids in which all 
particles except for one are coupled to the sites of a lattice via harmonic springs. The 
remaining particle is fixed in space and is called the carrier of the lattice. We fix this 
particle to the origin of the coordinate system. As we are interested in anisotropic particles, 
we will also restrict their rotational motion by a contribution U rot (6f) to the potential. We 
set this to be proportional to sin 2 ^, where 0$ is the angle between the axis of particle ri; 
and a reference axis m; (cf. Fig [2]). 
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Figure 2: Scheme of an Einstein Molecule for hard ellipsoids 



This kind of model is known as an Einstein Molecule (EM) 17] . (The reason for fixing 
one particle is the following: In the case of an Einstein Crystal (EC), center of mass motion 
of the entire system does not cost energy. Hence, for weak coupling one needs to shift all 
particle positions after every move to keep the center of mass positioned, as it was done 
e.g. in the work of Poison et al. [18( . In the case of the EM, the fixed carrier particle ensures 
non- divergency of the center of mass mean square displacement for a negligible harmonic 



potential (see also 



The interaction potential of the EM is 

^ref = J2' A * mnS ( r ' - T ^) 2 + Y! ^ Sin2 9i 

i i 

= X Ys' [(n-r ,,) 2 + sin 2 ^ 



(4) 



where r 0i j denote the position vectors of the lattice sites. The prime denotes that the sum 
runs over all particles except for the carrier. For simplicity we chose the spring constants 
of all lattice sites in the second line of Eq. H] as equal. In addition we set A trans = A rot = A. 
We use twice the short axis b as the unit of length and k^T as the unit of energy (except 
where stated otherwise). With this the unit of A is ksT / (2b) 2 . As we are only interested in 
the configurational part of phase space, the kinetic energies of the particles are disregarded. 
In order to evaluate the configurational part of the partition function of the Einstein 
Molecule, we assume that the maximum coupling constant A max is strong enough for 0j <C 1. 
So we obtain (cf. 
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(5) 



In case of the SM2-EM the same approach leads to 
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(6) 



The derivation of Eq. [6] is outlined in Appendix |A] The difference in free energy per particle 
between the FCC-EM and the SM2-EM is (In 2)/N, due to the presence of two types of lattice 
sites in the SM2 unit cell. This difference vanishes in the thermodynamic limit — > oo. 
Coming back to the calculation of F sys we rewrite the integral in Eq. [TJ (with ( = A/A max ) as 

-o 



AFo 



dX 



J2' sin 2 ^. 



(7) 



We evaluate eqs. [3] and [7] by the following expanded ensemble technique: We discretize 
the range of values for A. Then, besides translational and rotational moves, we perform a 
move in which the system passes from one model with a value A, into an adjacent one with Xj 
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and vice versa. In order to ensure good statistics when sampling the A-range, we introduce 
a set of weights ip m and sample the expanded ensemble given by the partition function 

M 

Z = ^Z m (A m )e^ m , (8) 

771=1 

where Z m (X m ) is the partition function of the model m with A = A m and ip m its weighting 
factor. The acceptance probability of a A-move is then given by 



Pi^j = min 



1 exp(^j - Uj) 
exp(^i - Ui) 



(9) 



such that for an adequate set of weights the system can be forced to visit the states of 
interest. One can then compute the free energy difference as 



Fj - Ft = - In 



Pi 



Here Pi and pj are the probabilities for the system to visit model i or model j, respectively, 
in the presence of the weights. 

The success of this procedure depends on finding appropriate weights. The weights are 
not known a priori, but they can be adjusted iteratively during the simulation, as has been 
introduced by Wang and Landau for the case of the density of states as a function of 
energy. We apply this idea to thermodynamic integration. Initially, we choose the weighting 
factors as ipi = Vz G {0, . . . , M}. Then simulations for each subensemble are carried out 
in which after each A-move the weight of the rejected model is increased by Aip = 1. This 
leads to an increase of both the possibility to visit the accepted model and the possibility 
to stay there. As in the original algorithm Aip is decreased by a factor a < 1, Ai/j — > a ■ Aip 
(here: a = 0.5), as soon as the difference between the probabilities becomes sufficiently 
small. Once the simulations are finished p m ~ p m -i (i n fact ln[p m /p m _i] was less than 
1(T 5 Nk B T after only 10 5 steps). 

Finally, we consider the computation of AFi = F off — F on , where the A-step does not 
change A but consists of switching on and off the hard-core potential. 

According to Eq. Amoves which switch off the potential or which lead to a state with no 
overlap are always accepted whereas moves of the form off — > on which yield a state with at 
least two overlaping particles are always rejected. For this case the coupling parameter was 
A = A max (i.e. the reference state), and hence the free energy difference between the states 
on and off was expected to be very small. 



Therefore we set the corresponding weights equal to and kept them fixed during the 
calculation. (This approach is validated by our results for Aiq, which were of the order of 
10~ 4 Nk B T.) 



III. RESULTS 



A. Hard spheres 

In order to test the algorithm before applying it to anisotropic particles, we first computed 
the free energy of hard spheres at various densities g = N/V and particle numbers N. Table 
[J summarizes our results. Fig. [3] shows the free energy per particle as a function oil/N for 
g = 1.04086. The dotted line is a fit to 

F , e2 , e3 



by which we extrapolate our results to infinte N (see also ref. 171]). 
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Figure 3: Hard spheres: Free energy per particle as a function of the inverse particle number. 
Symbols: MC data. Line: Fit according to eq. [11] 

Fig. H] shows a comparison of our results with Einstein Crystal and Einstein Molecule 
computations that did not use the Wang Landau algorithm. There is very good agreement, 
giving us confidence in the results for the hard ellipsoid case. 
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Figure 4: Hard sp 
methods ((a) ref. 



reres: Difference AF between our results and free energies computed by similar 



13] and (b) ref. 18]). There is very good agreement. 
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Table I: Results for the free energy of hard spheres 

N 

504 

768 

810 
1728 
1728 
1728 



F/N 
4.924(12) 
4.957(10) 
4.962(10) 
4.988(7) 
5.647(7) 
6.283(7) 



B. Hard ellipsoids 

Previous work showed that the angle of inclination of the SM2 unit cell, (3, is a very "soft" 
degree of freedom 12], i.e. the corresponding shear modulus is almost zero. (3 fluctuated 
strongly even at a pressure as high as P = 46 kBT/8a6 2 (for a/b = 3, where the nematic-solid 
coexistence pressure is P = 31 kBT/8a& 2 [1]). This unusual mechanical property is due to 
the fact that planes of equally oriented particles can slide across each other without much 
interaction, unless the system is forced to pack very densely. In order to quantify this effect, 
we computed free energies for various fixed values of (3. In the special case that the unit cell 
is invariant under reflections with respect to the bc-plane (see Fig. [T|), the configuration has 
the same symmetry as (but different unit cell parameters than) the close-packed structure 



constructed by Donev et al. 11] . In the following we refer to this structure as SM2^ cp \ 
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Fig. [5] shows the free energy for a/6 = 3 as a function of the density g (to simplify the 
comparison with other studies we use l/8a& 2 as the unit of density here instead of 1/86 3 ). 
The symbols are direct simulation results of free energies: triangles for SM2^ cp \ circles for 
SM2 with different values of f3 (see table [IT] for details) and a square for stretched fee. The 
solid lines are polynomial fits to the equation of state data from our previous work (ref. jl^). 
integrated over g and shifted by a constant to fit the free energy data. 

12i 1 1 1 1 1 1 1 — i 1 1 — i — ' — i — ' — i 
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Figure 5: Hard ellipsoids, a/b = 3: Free energy per particle as a function of density. Symbols are 
direct MC results for the free energy, lines are fits to MC data for the equation of state. 



Taking the errors into account there is no evidence for a difference in free energy be- 
tween the different angles of inclination (3 for the SM2 crystals. This supports our earlier 
observation that the angle of inclination is a soft degree of freedom [l^ . 

For decreasing density, the free energy difference between stretched fee and SM2 decreases 
and the lines intersect at g ~ 1.17, which is very close to the solid-nematic phase transition 
(g = 1.163 according to ref. [l|]). Our data therefore confirm that SM2 is more stable than 
fee at a/b = 3 and above g ~ 1.17. 

In Table 1111} we compare the free energies of SM2, SM2^ and fee function of 
aspect ratio, viz. for a/b = 1.55,2 and 3. As our input configurations were produced at 
a fixed pressure (P = 46 /cbT/8o6 2 ), systems of different aspect ratios and/or structure 
had different densities. In order to compare them, we calculated the Gibbs free energy per 
particle G/N by the Legendre transform of F/N with respect to the volume. Fig. [6] shows 
the Gibbs free energy per particle. (The lines are guides to the eye, only.) Again there is 
no difference (within the errorbars) between the SM2^ structure and the other values of 
(3. The superior stability of SM2 is confirmed for a/b > 2. At a/b = 1.55, however, we find 
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Table II: Hard Ellipsoids: Free energy per particle for a/b = 3 



Lattice type 



N 



0\ 



F sim /N 



FCC 
SM2 

SM2( c p) 



1.23390 
1.22545 
1.27111 
1.22815 
1.27352 



432 
432 
432 
432 
432 



(90) 
110.76 
115.91 
148.35 
147.97 



8.95(8) 
8.52(16) 
9.68(16) 
8.71(16) 
9.76(16) 



that fee is more stable, indicating a phase transition between a/b = 1.55 and a/b = 2.0. 

(This happens to be near a/b = a/3 , the lower boundary of aspect ratios for which 
prolate ellipsoids can form crystals with maximal packing fraction r\ = 0.770732 11]; but 
smaller aspect ratios near this value still produce higher-than-fee densities, so that we do 
not suspect a connection.) 
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Figure 6: Hard ellipsoids: Gibbs free energy per particle vs. aspect ratio a/b at P = 46/csT '/Sab 2 . 
Lines to guide the eye. At a/b > 2, SM2 is more stable, while fee is more stable at a/b = 1.55, 
implying a solid-solid phase transition in between. 

In Fig. [7] we show an updated phase diagram of hard ellipsoids of revolution. It includes 
part of the results of Frenkel and Mulder [l], and their suggested phase boundaries and 
coexistence regions. We have inserted our state points (this work and [7]), and extended 
its high-density boundary to the maximum densities found by Donev et al. [ll], hence 
including all densities possible in SM2 (recall that SM2 at maximum packing coincides with 
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Table III: Hard ellipsoids: Free energy and Gibbs free energy per particle (P = 46 IzbT /Rab 2 ) 
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the packings of Donev et al. ). As stated above, our data imply a phase transition between 
SM2 and fee near a/b = \^3. 

In hashes we indicate a possible location of the coexistence region, according to the 
following argument: For spheres (a/b = 1) and maximum packing (g = \/2) the density 
differences among plastic solid, fee and SM2 vanish, so the coexistence regions among these 
phases should join and vanish in width either at this point, or before this point is reached. 
(We can not make statements yet about the details of the approach to the sphere limit, 
therefore we have not drawn anything) . For a/b > 1, and above g = \/2 (the dotted 
line), only SM2 exists, so the SM2-fcc coexistence region must lie below g = The 
packing efficiency of SM2 and the resulting entropic advantage should favor SM2 even below 
g = v2, and the stronger this advantage, the lower the transition density - hence the 
downward slope of the coexistence region with increasing a/b. The packing advantage also 
dictates the increase in width of the region, since SM2 is accordingly higher in density at 
a given pressure. Finally, the coexistence region should pass between our state points of 
fee at a/b = 1.55, and SM2 at a/b = 2.00. The width there we estimate from the density 
difference at a/b = 1.55 and a/b = 2.00 at pressure P = 46 kBT/8ab 2 (Table [TTT]) . 
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Figure 7: (color online). Updated phase diagram of hard ellipsoids of revolution. It includes part 
of the results of Frenkel and Mulder [l[ (open symbols), and their suggested phase boundaries and 
coexistence regions. The data points at a/b = 1 are taken from [20]. We have inserted our state 
points (this work and [7]; filled symbols), and extended its high-density boundary to the maximum 



densities found by Donev et al. hence including all densities possible in SM2 (recall that SM2 
at maximum packing coincides with the packings of Donev et al. ) . In hashes we indicate a possible 
location of the coexistence region between fee and SM2 (see text for details). 



IV. CONCLUSION 



In summary, we have studied crystalline phases of hard ellipsoids considering their relative 
stability. We calculated the absolute free energies as functions of the particle density g 
and the aspect ratio a/b by use of a thermodynamic integration technique with an Einstein 
Molecule as the reference state. The integration path was sampled by an expanded ensemble 
method in which the weights were adjusted by the Wang-Landau algorithm. After checking 
our simulations for reliability considering the test case of hard spheres, we applied our 
methods to ellipsoids. At pressure P = 46 ksT/8ab 2 we found no difference in the free 
energies of SM2 crystals with different angles of inclination (3 . Furthermore our results 
show that the SM2 phase is more stable than the stretched-fee phase for densities g > 1.17 
(at a/b = 3) and for aspect ratios a/b > 2.0 (at P = 46 ksT /&ab 2 ). Hard ellipsoids exhibit 
a fcc-SM2 phase transition between a/b = 1.55 and a/b = 2.0. 
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Appendix A: CALCULATION OF F SM2 

First, without loss of generality, we label particle i — 1 as the carrier of the lattice. Then 
we write down the partition function of the SM2-EM using Eq. HJ 

Z SM2 = T (N) x I dj>i I ^ 



The two trivial integrations are due to our freedom of choosing ri as origin of the coordinate 
system and Q\ as some orientation in space. 

is a combinatorial factor: We consider a lattice G which consists of N particles with 
two different orientations (distinguished by the primes in Fig. [8]). We now divide G up into 
two sublattices G' and G" with respect to the particle types. On these sublattices there are 
(N/2)\ possibilities, respectively, to position the particles on their sites. To account for the 
presence of the carrier on one of the sublattices the associated factorial is (N/2 — 1)!. Hence 



The integral over the spatial coordinates can directly be carried out and leads to 
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Figure 8: Combinatorics for an Einstein Molecule of the SM2 type 



Here, because of the azimuthal symmetry of the problem, we already carried out the 
integration over and used the approximation sin# ss 6 which was motivated before. The 
two remaining integrals now are related to the sublattices G' and G" . Solving them, we find 
for the resulting partition function 
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